Open Case Studies : CO2 Emissions Over Time

Disclaimer: The purpose of the Open Case Studies project is to demonstrate the use of various data science methods, tools, and software in the context of messy, real-world data. A given case study does not cover all aspects of the research process, is not claiming to be the most appropriate way to analyze a given data set, and should not be used in the context of making policy decisions without external consultation from scientific experts.

Motivation

---

This case study explores how different countries have contributed to Carbon Dioxide (CO2) emissions over time and how CO2 emission rates may relate to increasing global temperatures and increased rates of natural disasters and storms. This report provides a basis for the motivation: https://www.epa.gov/report-environment/greenhouse-gases.

CO2 makes up the largest proportion of greenhouse gas emissions in the United States:

A variety of sources and sectors contribute to greenhouse gas emissions, with transportation contributing the most metric tons of CO2:

So why should we pay attention to greenhouse gases?

According to the US Environmental Protection Agency (EPA) Inventory of U.S. Greenhouse Gas Emissions and Sinks 2020 Report:

Greenhouse gases absorb infrared radiation, thereby trapping heat in the atmosphere and making the planet warmer. The most important greenhouse gases directly emitted by humans include carbon dioxide (CO2), methane (CH4), nitrous oxide (N2O), and several fluorine-containing halogenated substances. Although CO2, CH4, and N2O occur naturally in the atmosphere, human activities have changed their atmospheric concentrations. From the pre- industrial era (i.e., ending about 1750) to 2018, concentrations of these greenhouse gases have increased globally by 46, 165, and 23 percent, respectively (IPCC 2013; NOAA/ESRL 2019a, 2019b, 2019c).

• IPCC stands for the Intergovernmental Panel on Climate Change

There are many signs that our planet is experiencing warmer temperatures:

The connection between greenhouse gas levels and global temperatures and the influence of increased global temperatures on human health are motivated by these reports:

Melillo, J.M., T.C. Richmond, and G.W. Yohe (eds.). 2014. Climate change impacts in the United States: The third National Climate Assessment. U.S. Global Change Research Program.

2020. “Inventory of US Greenhouse Gas Emissions and Sinks: 1990–2018.” EPA 430-R-20-002, Tech. Rep. https://www.epa.gov/ghgemissions/inventory-us-greenhouse-gas-emissions-and-sinks.

The National Climate Assessment Report states that:

Heat-trapping gases already in the atmosphere have committed us to a hotter future with more climate-related impacts over the next few decades. The magnitude of climate change beyond the next few decades depends primarily on the amount of heat-trapping gases that human activities emit globally, now and in the future.

See here and here for more information.

Main Questions

---

Our main question:

1. How have global CO2 emission rates changed over time? In particular for the US, and how does the US compare to other countries?
2. Are US CO2 emissions, global temperatures, and US storm rates associated?

Learning Objectives

---

In this case study, we will explore CO2 emission data from around the world. We will also focus on the US specifically to evaluate patterns of temperatures and storm activity. This case study will particularly focus on visualizations of patterns over time. We will especially focus on using packages and functions from the Tidyverse, such as plotlyand gganimate. The tidyverse is a library of packages created by RStudio. While some students may be familiar with previous R programming packages, these packages make data science in R especially efficient.

---

We will begin by loading the packages that we will need:

library(here)
library(tidyverse)
library(readxl)
library(plotly)
library(gganimate)
library(RColorBrewer)

Package	Use
here	to easily load and save data
tidyverse	to wrangle the data and create ggplot2 plots
readxl	to import the Excel file data
plotyly	to make the visualizations
gganimate	to make the plots interactive
RColorBrewer	to have greater control over the color in our plots

The first time we use a function, we will use the :: to indicate which package we are using. Unless we have overlapping function names, this is not necessary, but we will include it here to be informative about where the functions we will use come from.

Context

---

Greenhouse gas emissions are due to both natural processes and anthropogenic (human-derived) activities.

These emissions are one of the contributing factors to rising global temperatures, which can have a great influence on public health as illustrated in the following image:

Gases in the atmosphere can contribute to climate change both directly and indirectly. Direct effects occur when the gas itself absorbs radiation. Indirect radiative forcing occurs when chemical transformations of the substance produce other greenhouse gases, when a gas influences the atmospheric lifetimes of other gases, and/or when a gas affects atmospheric processes that alter the radiative balance of the earth (e.g., affect cloud formation or albedo). The IPCC developed the Global Warming Potential (GWP) concept to compare the ability of a greenhouse gas to trap heat in the atmosphere relative to another gas. The GWP of a greenhouse gas is defined as the ratio of the accumulated radiative forcing within a specific time horizon caused by emitting 1 kilogram of the gas, relative to that of the reference gas CO2 (IPCC 2013). Therefore GWP-weighted emissions are provided in million metric tons of CO2 equivalent (MMT CO2 Eq.)

CO2 is actually the least capable of the greenhouse gases for trapping heat:

However, because CO2 is so much more abundant and stays in the atmosphere so much longer than other greenhouse gases, it has been the largest contributor to global warming.

See here for more details.

Furthermore, sizing CO2 levels also influence ocean acidity:

This makes it difficult for organisms to maintain their shells or skeletons that are made of calcium carbonate, thus making it more difficult for these organisms to survive and impacting their role in the ecosystem and food chain.

Furthermore, greenhouse gas emissions are believed to influence storm rates.

Indeed events with high levels of precipitation which can induce flooding and property damage are generally increasing around the country:

https://www.epa.gov/climate-indicators/cherry-blossoms

Limitations

---

There are some important considerations regarding this data analysis to keep in mind:

1. Limitation 1

2. Limitaiton 2

What are the data?

---

In this case study we will be using data related to CO2 emissions, as well as other data that may influence, be influenced or relate to CO2 emissions. Most of our data was obtained from Gapminder, which is a unique nonprofit that provides a variety of data for free.

In their words, Gapminder is…

Gapminder is an independent Swedish foundation with no political, religious or economic affiliations. Gapminder is a fact tank, not a think tank. Gapminder fights devastating misconceptions about global development. Gapminder produces free teaching resources making the world understandable based on reliable statistics. Gapminder promotes a fact-based worldview everyone can understand. Gapminder collaborates with universities, UN, public agencies and non-governmental organizations. All Gapminder activities are governed by the board. We do not award grants. Gapminder Foundation is registered at Stockholm County Administration Board. Our constitution can be found here.

Data	Source	Orginal Source	Description
CO2 emissions	Gapminder	Carbon Dioxid Information Analysis Center (CDIAC)	CO2 emissions in metric tonnes per person
GDP per capita, yearly growth	Gapminder	World Bank	Growth Domestic Product (which is a measure of economic status) per person
Energy use per person	Gapminder	World Bank	Use of primary energy before transformation to other end-use fules, by country starting in 1960
Crude Mortality Rate	World Bank	World Bank	Death rate per 1,000 people from 1960 to 2018 by country
US Natural Disasters	The National Oceanic and Atmospheric Administration (NOAA)	The National Oceanic and Atmospheric Administration (NOAA)	Data from 1980 to 2019 including: – Drought Count – Flooding Count – Freeze Count – Severe Storm Count – Tropical Cyclone Count – Wildfire Count – Winter Storm Count
Temperature	The National Oceanic and Atmospheric Administration (NOAA)	The National Oceanic and Atmospheric Administration (NOAA)	National yearly average temperature (in Fahrenheit) from 1980 to 2019

not quite right… almost… https://www.ncdc.noaa.gov/cag/national/time-series/110/tavg/12/5/1895-2020?base_prd=true&begbaseyear=1901&endbaseyear=2000

https://www.ncdc.noaa.gov/cag/national/time-series/110/pcp/12/12/1895-2020?base_prd=true&begbaseyear=1901&endbaseyear=2000

I think this is what michael used: https://www.ncdc.noaa.gov/cag/national/time-series/110/tavg/1/12/1895-2020?base_prd=true&begbaseyear=1901&endbaseyear=2000

No… not it either…

variable2 | Variable info
– more details
– more detials
Example: Content content | source | Description

Data Import

---

CO2_emissions <- read_excel(here("docs/yearly_co2_emissions_1000_tonnes.xlsx"))
gdp_growth <- read_xlsx(here("docs/gdp_per_capita_yearly_growth.xlsx"))
energy_use <- read_xlsx(here("docs/energy_use_per_person.xlsx"))
mortality <- read_xls(here("docs/API_SP.DYN.CDRT.IN_DS2_en_excel_v2_804384.xls"))

us_disaster <- read_csv(here("docs/time-series-US.csv"), skip = 2)
us_temperature <- read_csv(here("docs/temperature.txt"),skip=4,na="-99")

Data Exploration and Wrangling

---

Global Data

We will also use the %>% pipe which can be used to define the input for later sequential steps. This will make more sense when we have multiple sequential steps using the same data object. To use the pipe notation we need to install and load dplyr as well.

Yearly CO₂ Emissions, 1000 Metric Tonnes

head(CO2_emissions)

# A tibble: 6 x 265
  country `1751` `1752` `1753` `1754` `1755` `1756` `1757` `1758` `1759` `1760`
  <chr>    <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1 Afghan…     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
2 Albania     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
3 Algeria     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
4 Andorra     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
5 Angola      NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
6 Antigu…     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
# … with 254 more variables: `1761` <dbl>, `1762` <dbl>, `1763` <dbl>,
#   `1764` <dbl>, `1765` <dbl>, `1766` <dbl>, `1767` <dbl>, `1768` <dbl>,
#   `1769` <dbl>, `1770` <dbl>, `1771` <dbl>, `1772` <dbl>, `1773` <dbl>,
#   `1774` <dbl>, `1775` <dbl>, `1776` <dbl>, `1777` <dbl>, `1778` <dbl>,
#   `1779` <dbl>, `1780` <dbl>, `1781` <dbl>, `1782` <dbl>, `1783` <dbl>,
#   `1784` <dbl>, `1785` <dbl>, `1786` <dbl>, `1787` <dbl>, `1788` <dbl>,
#   `1789` <dbl>, `1790` <dbl>, `1791` <dbl>, `1792` <dbl>, `1793` <dbl>,
#   `1794` <dbl>, `1795` <dbl>, `1796` <dbl>, `1797` <dbl>, `1798` <dbl>,
#   `1799` <dbl>, `1800` <dbl>, `1801` <dbl>, `1802` <dbl>, `1803` <dbl>,
#   `1804` <dbl>, `1805` <dbl>, `1806` <dbl>, `1807` <dbl>, `1808` <dbl>,
#   `1809` <dbl>, `1810` <dbl>, `1811` <dbl>, `1812` <dbl>, `1813` <dbl>,
#   `1814` <dbl>, `1815` <dbl>, `1816` <dbl>, `1817` <dbl>, `1818` <dbl>,
#   `1819` <dbl>, `1820` <dbl>, `1821` <dbl>, `1822` <dbl>, `1823` <dbl>,
#   `1824` <dbl>, `1825` <dbl>, `1826` <dbl>, `1827` <dbl>, `1828` <dbl>,
#   `1829` <dbl>, `1830` <dbl>, `1831` <dbl>, `1832` <dbl>, `1833` <dbl>,
#   `1834` <dbl>, `1835` <dbl>, `1836` <dbl>, `1837` <dbl>, `1838` <dbl>,
#   `1839` <dbl>, `1840` <dbl>, `1841` <dbl>, `1842` <dbl>, `1843` <dbl>,
#   `1844` <dbl>, `1845` <dbl>, `1846` <dbl>, `1847` <dbl>, `1848` <dbl>,
#   `1849` <dbl>, `1850` <dbl>, `1851` <dbl>, `1852` <dbl>, `1853` <dbl>,
#   `1854` <dbl>, `1855` <dbl>, `1856` <dbl>, `1857` <dbl>, `1858` <dbl>,
#   `1859` <dbl>, `1860` <dbl>, …

CO2_emissions <- CO2_emissions %>%
  gather(key = Year, value = Emissions, -country) %>%
  rename(Country=country) %>%
  mutate(Emissions = Emissions/1000) %>%
  rename(`CO2 Emissions (Mg)`=Emissions)

sapply(CO2_emissions, class)

           Country               Year CO2 Emissions (Mg) 
       "character"        "character"          "numeric" 

CO2_emissions$Year <- as.numeric(CO2_emissions$Year)

summary(CO2_emissions$Year)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1751    1817    1882    1882    1948    2014 

Yearly Growth in GDP per Capita

Michael Need to determine units for this

head(gdp_growth)

# A tibble: 6 x 220
  country   `1801`   `1802`   `1803`   `1804`   `1805`   `1806`   `1807`
  <chr>      <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
1 Afghan… NA       NA       NA       NA       NA       NA       NA      
2 Albania  0.104    0.104    0.104    0.104    0.104    0.104    0.104  
3 Algeria -0.00247 -0.00247 -0.00247 -0.00247 -0.00247 -0.00247 -0.00247
4 Andorra  0.166    0.166    0.166    0.166    0.166    0.166    0.166  
5 Angola   0.425    0.425    0.425    0.425    0.425    0.425    0.425  
6 Antigu… NA       NA       NA       NA       NA       NA       NA      
# … with 212 more variables: `1808` <dbl>, `1809` <dbl>, `1810` <dbl>,
#   `1811` <dbl>, `1812` <dbl>, `1813` <dbl>, `1814` <dbl>, `1815` <dbl>,
#   `1816` <dbl>, `1817` <dbl>, `1818` <dbl>, `1819` <dbl>, `1820` <dbl>,
#   `1821` <dbl>, `1822` <dbl>, `1823` <dbl>, `1824` <dbl>, `1825` <dbl>,
#   `1826` <dbl>, `1827` <dbl>, `1828` <dbl>, `1829` <dbl>, `1830` <dbl>,
#   `1831` <dbl>, `1832` <dbl>, `1833` <dbl>, `1834` <dbl>, `1835` <dbl>,
#   `1836` <dbl>, `1837` <dbl>, `1838` <dbl>, `1839` <dbl>, `1840` <dbl>,
#   `1841` <dbl>, `1842` <dbl>, `1843` <dbl>, `1844` <dbl>, `1845` <dbl>,
#   `1846` <dbl>, `1847` <dbl>, `1848` <dbl>, `1849` <dbl>, `1850` <dbl>,
#   `1851` <dbl>, `1852` <dbl>, `1853` <dbl>, `1854` <dbl>, `1855` <dbl>,
#   `1856` <dbl>, `1857` <dbl>, `1858` <dbl>, `1859` <dbl>, `1860` <dbl>,
#   `1861` <dbl>, `1862` <dbl>, `1863` <dbl>, `1864` <dbl>, `1865` <dbl>,
#   `1866` <dbl>, `1867` <dbl>, `1868` <dbl>, `1869` <dbl>, `1870` <dbl>,
#   `1871` <dbl>, `1872` <dbl>, `1873` <dbl>, `1874` <dbl>, `1875` <dbl>,
#   `1876` <dbl>, `1877` <dbl>, `1878` <dbl>, `1879` <dbl>, `1880` <dbl>,
#   `1881` <dbl>, `1882` <dbl>, `1883` <dbl>, `1884` <dbl>, `1885` <dbl>,
#   `1886` <dbl>, `1887` <dbl>, `1888` <dbl>, `1889` <dbl>, `1890` <dbl>,
#   `1891` <dbl>, `1892` <dbl>, `1893` <dbl>, `1894` <dbl>, `1895` <dbl>,
#   `1896` <dbl>, `1897` <dbl>, `1898` <dbl>, `1899` <dbl>, `1900` <dbl>,
#   `1901` <dbl>, `1902` <dbl>, `1903` <dbl>, `1904` <dbl>, `1905` <dbl>,
#   `1906` <dbl>, `1907` <dbl>, …

gdp_growth <- gdp_growth %>%
  gather(key = Year, value = gdp_growth, -country) %>%
  rename(Country=country) %>%
  rename(`GDP Growth/Capita (%)` = gdp_growth) %>%
  filter(complete.cases(.))

sapply(gdp_growth, class)

              Country                  Year GDP Growth/Capita (%) 
          "character"           "character"             "numeric" 

gdp_growth$Year <- as.numeric(gdp_growth$Year)

summary(gdp_growth$Year)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1801    1863    1915    1915    1967    2019 

Energy Use per Person

head(energy_use)

# A tibble: 6 x 57
  country `1960` `1961` `1962` `1963` `1964` `1965` `1966` `1967` `1968` `1969`
  <chr>    <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1 Albania     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
2 Algeria     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
3 Angola      NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
4 Antigu…     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
5 Argent…     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
6 Armenia     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
# … with 46 more variables: `1970` <dbl>, `1971` <dbl>, `1972` <dbl>,
#   `1973` <dbl>, `1974` <dbl>, `1975` <dbl>, `1976` <dbl>, `1977` <dbl>,
#   `1978` <dbl>, `1979` <dbl>, `1980` <dbl>, `1981` <dbl>, `1982` <dbl>,
#   `1983` <dbl>, `1984` <dbl>, `1985` <dbl>, `1986` <dbl>, `1987` <dbl>,
#   `1988` <dbl>, `1989` <dbl>, `1990` <dbl>, `1991` <dbl>, `1992` <dbl>,
#   `1993` <dbl>, `1994` <dbl>, `1995` <dbl>, `1996` <dbl>, `1997` <dbl>,
#   `1998` <dbl>, `1999` <dbl>, `2000` <dbl>, `2001` <dbl>, `2002` <dbl>,
#   `2003` <dbl>, `2004` <dbl>, `2005` <dbl>, `2006` <dbl>, `2007` <dbl>,
#   `2008` <dbl>, `2009` <dbl>, `2010` <dbl>, `2011` <dbl>, `2012` <dbl>,
#   `2013` <dbl>, `2014` <dbl>, `2015` <dbl>

energy_use <- energy_use %>%
  gather(key = Year, value = energy_use, -country) %>%
  rename(Country=country) %>%
  rename(`Energy Use (kg, oil-eq./capita)` = energy_use)

sapply(energy_use, class)

                        Country                            Year 
                    "character"                     "character" 
Energy Use (kg, oil-eq./capita) 
                      "numeric" 

energy_use$Year <- as.numeric(energy_use$Year)

summary(energy_use$Year)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1960    1974    1988    1988    2001    2015 

Crude Mortality Rate

head(mortality)

# A tibble: 6 x 64
  `Data Source` `World Developm… ...3  ...4  ...5  ...6  ...7  ...8  ...9  ...10
  <chr>         <chr>            <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Last Updated… 43819            <NA>  <NA>  <NA>  <NA>  <NA>  <NA>  <NA>  <NA> 
2 <NA>          <NA>             <NA>  <NA>  <NA>  <NA>  <NA>  <NA>  <NA>  <NA> 
3 Country Name  Country Code     Indi… Indi… 1960  1961  1962  1963  1964  1965 
4 Aruba         ABW              Deat… SP.D… 6.38… 6.24… 6.11… 6.01… 5.91… 5.83…
5 Afghanistan   AFG              Deat… SP.D… 32.2… 31.6… 31.0… 30.5… 30.0… 29.5…
6 Angola        AGO              Deat… SP.D… 27.0… 26.8… 26.6… 26.4… 26.1… 25.9…
# … with 54 more variables: ...11 <chr>, ...12 <chr>, ...13 <chr>, ...14 <chr>,
#   ...15 <chr>, ...16 <chr>, ...17 <chr>, ...18 <chr>, ...19 <chr>,
#   ...20 <chr>, ...21 <chr>, ...22 <chr>, ...23 <chr>, ...24 <chr>,
#   ...25 <chr>, ...26 <chr>, ...27 <chr>, ...28 <chr>, ...29 <chr>,
#   ...30 <chr>, ...31 <chr>, ...32 <chr>, ...33 <chr>, ...34 <chr>,
#   ...35 <chr>, ...36 <chr>, ...37 <chr>, ...38 <chr>, ...39 <chr>,
#   ...40 <chr>, ...41 <chr>, ...42 <chr>, ...43 <chr>, ...44 <chr>,
#   ...45 <chr>, ...46 <chr>, ...47 <chr>, ...48 <chr>, ...49 <chr>,
#   ...50 <chr>, ...51 <chr>, ...52 <chr>, ...53 <chr>, ...54 <chr>,
#   ...55 <chr>, ...56 <chr>, ...57 <chr>, ...58 <chr>, ...59 <chr>,
#   ...60 <chr>, ...61 <chr>, ...62 <chr>, ...63 <chr>, ...64 <chr>

colnames(mortality) <- mortality[3,]
mortality <- mortality[-c(1:3),]

mortality <- mortality %>%
  select(-`Country Code`,
         -`Indicator Name`,
         -`Indicator Code`) %>%
  rename(Country="Country Name") %>%
  gather(key = Year, value = `Deaths/1000 People`, -Country)

sapply(mortality, class)

           Country               Year Deaths/1000 People 
       "character"        "character"        "character" 

mortality$Year <- as.numeric(mortality$Year)
mortality$`Deaths/1000 People` <- as.numeric(mortality$`Deaths/1000 People`)

summary(mortality$Year)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1960    1975    1990    1990    2004    2019 

US-specific Data

Disasters

us_disaster <- us_disaster %>%
  dplyr::select(Year,
                `Drought Count`,
                `Flooding Count`,
                `Freeze Count`,
                `Severe Storm Count`,
                `Tropical Cyclone Count`,
                `Wildfire Count`,
                `Winter Storm Count`)

sapply(us_disaster, class)

                  Year          Drought Count         Flooding Count 
             "numeric"              "numeric"              "numeric" 
          Freeze Count     Severe Storm Count Tropical Cyclone Count 
             "numeric"              "numeric"              "numeric" 
        Wildfire Count     Winter Storm Count 
             "numeric"              "numeric" 

us_disaster <- us_disaster %>%
  mutate(`Disasters` = rowSums(.[2:dim(.)[2]])) %>%
  dplyr::select(Year,
                `Disasters`) %>%
  mutate(Country="United States") %>%
  gather(key = Indicator, value = Value, -Country, -Year) %>%
  mutate(Region="United States",
         Type="US-specific")

Temperature

us_temperature <- read_csv(here("docs/temperature.txt"),skip=4,na="-99")

sapply(us_temperature, class)

     Date     Value   Anomaly 
"numeric" "numeric" "numeric" 

us_temperature <- us_temperature %>%
  dplyr::select(-Anomaly) %>%
  mutate(Date = substr(Date, start = 1, stop = 4)) %>%
  mutate(Country="United States",
         Year = Date,
         Indicator="Temperature (Fahrenheit)",
         Region="United States",
         Type="US-specific") %>%
  dplyr::select(-Date, )

Analysis data frame

colnames(CO2_emissions)

[1] "Country"            "Year"               "CO2 Emissions (Mg)"

colnames(gdp_growth)

[1] "Country"               "Year"                  "GDP Growth/Capita (%)"

colnames(energy_use)

[1] "Country"                         "Year"                           
[3] "Energy Use (kg, oil-eq./capita)"

colnames(mortality)

[1] "Country"            "Year"               "Deaths/1000 People"

df_wide <- CO2_emissions %>%
  full_join(gdp_growth, by=c("Country", "Year")) %>%
  full_join(energy_use, by=c("Country", "Year")) %>%
  full_join(mortality, by=c("Country", "Year"))

df_long <- df_wide %>%
  gather(key=Indicator,value=Value,-Country, -Year) %>%
  mutate(Region=case_when(Country=="United States" ~ "United States",
                          Country!="United States" ~ "Rest of the World"),
         Type="Global")

setequal(sapply(df_long, class),sapply(us_disaster, class))

[1] TRUE

setequal(sapply(df_long, class),sapply(us_temperature, class))

[1] TRUE

df_long <- df_long %>%
  rbind(us_disaster) %>%
  rbind(us_temperature)

sapply(df_long,class)

    Country        Year   Indicator       Value      Region        Type 
"character" "character" "character"   "numeric" "character" "character" 

df_long$Country <- as.factor(df_long$Country)
df_long$Year <- as.numeric(df_long$Year)

sapply(df_long,class)

    Country        Year   Indicator       Value      Region        Type 
   "factor"   "numeric" "character"   "numeric" "character" "character" 

df_long <- df_long %>%
  filter(complete.cases(.)) %>%
  arrange(Country)

Data Visualization

---

CO2 plots

ggplot(CO2_emissions, aes(x=Year, y=`CO2 Emissions (Mg)`, group=Country)) +
  geom_line(alpha=0.2) + 
  labs(title = expression("Country CO"[2]*" Emissions per Year , 1751-2014"),
       caption = expression("Limited to reporting countries")) + 
  ylab("Emissions (1M Metric Tonnes)")

CO2_world<-CO2_emissions %>%
  group_by(Year) %>%
  summarise(`CO2 Emissions (Mg)` = sum(`CO2 Emissions (Mg)`, na.rm = TRUE)) %>%
  ggplot(aes(x=Year, y=`CO2 Emissions (Mg)`)) +
  geom_line() + 
  labs(title = expression("World CO"[2]*" Emissions per Year , 1751-2014"),
       caption = expression("Limited to reporting countries")) + 
  ylab("Emissions (1M Metric Metric Tonnes)")
CO2_world

Yearly Growth in GDP per Capita

ggplot(gdp_growth, aes(x=Year, y=`GDP Growth/Capita (%)`, group=Country)) +
  geom_line(alpha=0.2) + 
  labs(title = expression("Country GDP Growth per Capita per Year (Annual %), 1801-2019"),
       caption = expression("Limited to reporting countries")) + 
  ylab("GDP Growth per Capita (Annual %)")

gdp_growth %>%
  group_by(Year) %>%
  summarise(`GDP Growth/Capita (%)` = mean(`GDP Growth/Capita (%)`, na.rm = TRUE)) %>%
  ggplot(aes(x=Year, y=`GDP Growth/Capita (%)`)) +
  geom_line() + 
  labs(title = expression("Mean Country GDP Growth per Capita per Year (Annual %), 1801-2019"),
       caption = expression("Limited to reporting countries")) + 
  ylab("GDP Growth per Capita (Annual %)")

Energy Use per Person

ggplot(energy_use, aes(x=Year, y=`Energy Use (kg, oil-eq./capita)`, group=Country)) +
  geom_line(alpha=0.2) + 
  labs(title = expression("Country Energy Use (kg of Oil Equivalent per Capita), 1960-2015",),
       caption = expression("Limited to reporting countries")) + 
  ylab("Energy Use (kg of Oil Equivalent per Capita)")

energy_use %>%
  group_by(Year) %>%
  summarise(`Energy Use (kg, oil-eq./capita)` = sum(`Energy Use (kg, oil-eq./capita)`, na.rm = TRUE)) %>%
  ggplot(aes(x=Year, y=`Energy Use (kg, oil-eq./capita)`)) +
  geom_line() + 
  labs(title = expression("Worldwide Energy Use (kg of Oil Equivalent per Capita), 1960-2015"),
       caption = expression("Limited to reporting countries")) + 
  ylab("Energy Use (kg of Oil Equivalent per Capita)")

Crude Mortality Rate

ggplot(mortality, aes(x=Year, y=`Deaths/1000 People`, group=Country)) +
  geom_line(alpha=0.2) + 
  labs(title = expression("Country Crude Mortality Rate (per 1000 Persons), 1960-2019",),
       caption = expression("Limited to reporting countries")) + 
  ylab("Crude Mortality Rate (per 1000 Persons)")

mortality %>%
  group_by(Year) %>%
  summarise(`Deaths/1000 People` = mean(`Deaths/1000 People`, na.rm = TRUE)) %>%
  ggplot(aes(x=Year, y=`Deaths/1000 People`)) +
  geom_line() + 
  labs(title = expression("Mean Country Crude Mortality Rate (per 1000 Persons), 1960-2019",),
       caption = expression("Limited to reporting countries")) + 
  ylab("Crude Mortality Rate (per 1000 Persons)")

US-specific Data

Disasters

ggplot(df_long, aes(x=Year, y=Value, group=Country)) +
  geom_line(alpha=0.2) + 
  facet_grid(Indicator~., scales = "free_y") +
  ylab("Indicator Value") + 
  labs(title="Distribution of Indicators by Year and Value")

Subsetting the data

df_long %>%
  filter(Type=="Global") %>%
  group_by(Year,Indicator) %>%
  tally() %>%
  ggplot(aes(x=Year, y=n, color=Indicator)) +
  geom_vline(xintercept = 1980, linetype=2, color="black") +
  geom_vline(xintercept = 2010, linetype=2, color="black") +
  geom_line() +
  labs(title = "Countries with Complete Data per Year",
       subtitle = "Global Data") + 
  ylab("Countries") + 
  scale_x_continuous(breaks = seq(1750,2020,by=10),
                     labels = seq(1750,2020,by=10)) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank(),
        legend.position = "bottom")

df_long %>%
  filter(Region=="United States") %>%
  group_by(Indicator) %>%
  summarise(Start=min(Year), End=max(Year)) %>%
  ggplot(aes(x=Indicator, y=End)) +
  geom_hline(yintercept = 1980, linetype=2, color="black") +
  geom_hline(yintercept = 2010, linetype=2, color="black") +
  geom_segment(aes(x=Indicator,
                   xend=Indicator,
                   yend=End,
                   y=Start)) +
  geom_point(aes(x=Indicator, y=Start), shape=16, color="black") +
  geom_point(aes(x=Indicator, y=End), shape=21, fill="white", color="black") + 
  coord_flip() +
  labs(title = "Complete Data per Year",
       subtitle = "US-specific Data") + 
  ylab("Countries") + 
  scale_y_continuous(breaks = seq(1750,2020,by=10),
                     labels = seq(1750,2020,by=10)) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90),
        axis.title = element_blank())

animation_1 <- df_long %>%
  filter(Type=="Global") %>%
  filter(Indicator=="Deaths/1000 People") %>%
  filter(Year>=1980) %>%
  filter(Year<=2010) %>%
  ggplot(aes(x=Year, y=Value, group=Country, color=Region, size=Region,alpha=Region)) +
  geom_point() +
  scale_color_manual(values = c("Red","Black")) +
  scale_alpha_manual(values = c(0.1, 1)) +
  scale_size_manual(values = c(0.25, 2)) +
  labs(title="Distribution of Indicators by Year and Value, 1980-2010") +
  theme(axis.text.x = element_text(angle = 90)) + 
  theme_classic() +
  ylab("Crude Mortality Rate") +
  transition_time(as.integer(Year)) +
  shadow_wake(wake_length = 1, alpha = FALSE)

animate(animation_1, fps = 10, duration = 5)

animation_2 <- df_long %>%
  filter(Type=="Global") %>%
  filter(Indicator=="Energy Use (kg, oil-eq./capita)") %>%
  filter(Year>=1980) %>%
  filter(Year<=2010) %>%
  ggplot(aes(x=Year, y=Value, group=Country, color=Region, size=Region, alpha=Region)) +
  geom_point() +
  scale_color_manual(values = c("Red","Black")) +
  scale_alpha_manual(values = c(0.1, 1)) +
  scale_size_manual(values = c(0.25, 2)) +
  labs(title="Distribution of Indicators by Year and Value, 1980-2010") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90)) + 
  ylab("Energy Use per Capita") +
  transition_time(as.integer(Year)) +
  shadow_wake(wake_length = 1, alpha = FALSE)

animate(animation_2, fps = 10, duration = 5)

animation_3 <- df_long %>%
  filter(Type=="Global") %>%
  filter(Indicator=="GDP Growth/Capita (%)") %>%
  filter(Year>=1980) %>%
  filter(Year<=2010) %>%
  ggplot(aes(x=Year, y=Value, group=Country, color=Region, size=Region, alpha=Region)) +
  geom_point() +
  scale_color_manual(values = c("Red","Black")) +
  scale_alpha_manual(values = c(0.1, 1)) +
  scale_size_manual(values = c(0.25, 2)) +
  labs(title="Distribution of Indicators by Year and Value, 1980-2010") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90)) +
  ylab("GDP Growth per Capita (%)") +
  transition_time(as.integer(Year)) +
  shadow_wake(wake_length = 1, alpha = FALSE)

animate(animation_3, fps = 10, duration = 5)

animation_4 <- df_long %>%
  filter(Type=="Global") %>%
  filter(Indicator=="CO2 Emissions (Mg)") %>%
  filter(Year>=1980) %>%
  filter(Year<=2010) %>%
  ggplot(aes(x=Year, y=Value, group=Country, color=Region, size=Region, alpha=Region)) +
  geom_point() +
  scale_color_manual(values = c("Red","Black")) +
  scale_alpha_manual(values = c(0.1, 1)) +
  scale_size_manual(values = c(0.25, 2)) +
  labs(title="Distribution of Indicators by Year and Value, 1980-2010") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90)) +
  ylab("CO2 Emissions (Mg)") +
  transition_time(as.integer(Year)) +
  shadow_wake(wake_length = 1, alpha = FALSE)

animate(animation_4, fps = 10, duration = 5)

Top10<-df_long %>%
  filter(Type=="Global") %>%
  filter(Indicator=="CO2 Emissions (Mg)") %>%
  filter(Year>=1900) %>%
  filter(Year<=2010) %>%
  group_by(Country) %>%
  mutate(max_val = max(Value)) %>%
  ungroup() %>%
  mutate(rank=dense_rank(-max_val)) %>%
  filter(rank<=10) %>%
  ggplot(aes(x=Year, y=fct_reorder(Country, Value, max))) +
  geom_tile(color="transparent", aes(fill=log(Value))) +
  scale_fill_gradientn(colors = c("yellow","red","black")) +
  scale_x_continuous(breaks = seq(1900,2010,by=5),
                     labels = seq(1900,2010,by=5)) + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90),
        axis.title = element_blank(),
        legend.position = "bottom") + 
  labs(title = "Top 10 Emissions-producing Countries in 2010 (1900-2010)",
       subtitle = "Ordered by Emissions Produced in 2010",
       fill = "Ln(CO2 Emissions (Mg))")

Top10

US-specific

df_long_us <- df_long %>%
  filter(Country=="United States")

# Approximated derivative function)
df_long_us <- df_long_us %>%
  filter(Year>=1900,
         Year<=2010) %>%
  group_by(Indicator) %>%
  mutate("Change (%)"=((Value/lag(Value))*100)-100,
         Mean=mean(Value),
         Anomaly=Value-Mean) %>%
  ungroup() %>%
  mutate(Anomaly_color=ifelse(Anomaly>0,"Positive",
                              ifelse(Anomaly<0,"Negative","Zero")),
         Anomaly_color=factor(Anomaly_color, levels = c("Positive",
                                                        "Negative",
                                                        "Zero"),
                              ordered = TRUE))

US_Indicators<-df_long_us %>%
  filter(Year>=1980) %>%
  filter(Year<=2010) %>%
  ggplot(aes(x=Year, y=Value)) + 
  geom_line() + 
  facet_wrap(Indicator~., ncol=2, nrow=3, scales = "free_y") + 
  scale_x_continuous(breaks = seq(1980,2010,by=5),
                     labels = seq(1980,2010,by=5)) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90),
        axis.title.y = element_blank()) + 
  labs(title = "US-specific Indicators (1980-2010)")
US_Indicators

df_long_us %>%
  filter(Year>=1980) %>%
  filter(Year<=2010) %>%
  ggplot(aes(x=Year, y=`Change (%)`, color=Indicator, fill="transparent")) + 
  geom_hline(yintercept=0.8, linetype=2) +
  geom_hline(yintercept=1.2, linetype=2) +
  geom_hline(yintercept = 1, linetype=3) +
  geom_line(size=0.5) + 
  scale_x_continuous(breaks = seq(1980,2010,by=5),
                     labels = seq(1980,2010,by=5)) +
  scale_y_continuous(breaks = seq(-500,1250, by=250),
                     labels = seq(-500,1250, by=250)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x  = element_blank(),
        legend.position = "bottom",
        legend.direction = "horizontal") + 
  labs(title = "US-specific Indicators (1980-2010)",
       subtitle = "Change (%) Lines")

df_long_us %>%
  filter(Year>=1980) %>%
  filter(Year<=2010) %>%
  filter(Indicator=="Temperature (Fahrenheit)"|
           Indicator=="Energy Use (kg, oil-eq./capita)"|
           Indicator=="CO2 Emissions (Mg)") %>%
  ggplot(aes(x=Year, y=`Change (%)`, color=Indicator)) + 
  geom_hline(yintercept=0.8, linetype=2) +
  geom_hline(yintercept=1.2, linetype=2) +
  geom_hline(yintercept = 1, linetype=3) +
  geom_line(size=1) +
  scale_x_continuous(breaks = seq(1980,2010,by=5),
                     labels = seq(1980,2010,by=5)) +
  scale_y_continuous(breaks = seq(-10,10, by=1),
                     labels = seq(-10,10, by=1),
                     limits = c(-10,10)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x  = element_blank(),
        legend.position = "bottom",
        legend.direction = "horizontal") + 
  labs(title = "Emissions, Energy Use, and Temperature (1980-2010)",
       subtitle = "Change (%) Lines")

df_long_us %>%
  filter(Year>=1980) %>%
  filter(Year<=2010) %>%
  filter(Indicator=="Temperature (Fahrenheit)"|
           Indicator=="Energy Use (kg, oil-eq./capita)"|
           Indicator=="CO2 Emissions (Mg)") %>%
  ggplot(aes(x=Year, y=`Change (%)`, color=Indicator)) + 
  geom_hline(yintercept=0.8, linetype=2) +
  geom_hline(yintercept=1.2, linetype=2) +
  geom_hline(yintercept = 1, linetype=3) +
  geom_smooth(size=1, alpha=0.1, aes(fill=Indicator), se=FALSE) +
  scale_x_continuous(breaks = seq(1980,2010,by=5),
                     labels = seq(1980,2010,by=5)) +
  scale_y_continuous(breaks = seq(-10,10, by=1),
                     labels = seq(-10,10, by=1),
                     limits = c(-10,10)) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x  = element_blank(),
        legend.position = "bottom",
        legend.direction = "horizontal") + 
  labs(title = "US Emissions, Energy Use, and Temperatures (1980-2010)",
       subtitle = "Smoothed Change (%) Lines")

df_long_us %>%
  filter(Year>=1980) %>%
  filter(Year<=2010) %>%
  filter(Indicator=="CO2 Emissions (Mg)"|
           Indicator=="Temperature (Fahrenheit)") %>%
  ggplot(aes(x=Year, y=Value)) +
  geom_line() + 
  scale_x_continuous(breaks = seq(1980,2010,by=5),
                     labels = seq(1980,2010,by=5)) + 
  facet_wrap(Indicator~., scales = "free_y", ncol=1) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90),
        axis.title = element_blank()) + 
  labs(title="US Emissions and Temperatures (1980-2010)")

df_long_us %>%
  filter(Indicator=="CO2 Emissions (Mg)"|
           Indicator=="Temperature (Fahrenheit)") %>%
  ggplot(aes(x=Year, y=Value)) +
  geom_vline(xintercept = 1980, linetype=2, color="black") + 
  geom_vline(xintercept = 2010, linetype=2, color="black") +
  geom_segment(aes(x=Year, y=Value, xend=Year, yend=Mean,color=Anomaly_color), size=1.25) +
  scale_color_manual(values = c("red","blue","gray")) + 
  geom_hline(aes(yintercept=Mean), linetype=1, color="black") +
  scale_x_continuous(breaks = seq(1900,2010,by=5),
                     labels = seq(1900,2010,by=5)) +
  facet_wrap(Indicator~., scales = "free_y", ncol=1) + 
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90),
        axis.title = element_blank(),
        legend.position = "none")  +
  labs(title = "US Emissions and Temperatures (1900-2010)",
       subtitle = "Indicator Mean Represented by Solid Black Line")

df_long_us %>%
  filter(Indicator=="Temperature (Fahrenheit)"|
           Indicator=="CO2 Emissions (Mg)") %>%
  ggplot(aes(x=Year, y=`Change (%)`)) +
  annotate("rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=Inf, alpha=0.25, fill="green") +
  annotate("rect", xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=0, alpha=0.25, fill="red") +
  geom_hline(yintercept=0, linetype=1) + 
  geom_segment(aes(x=Year, y=`Change (%)`, xend=Year, yend=0), size=1.25) +
  facet_wrap(Indicator~., scales = "free_y", ncol=1) +
  scale_x_continuous(breaks = seq(1900,2010,by=5),
                     labels = seq(1900,2010,by=5)) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90),
        axis.title.x  = element_blank(),
        legend.position = "bottom",
        legend.direction = "horizontal") + 
  labs(title = "US Emissions, Energy Use, and Temperatures (1980-2010)",
       subtitle = "Change (%) Lines")

Main plot

library(patchwork)

CO2_world + Top10 + US_Indicators +
  plot_layout(widths = c(1, 2), heights = unit(c(2, 5), c('cm', 'null')))

png(here::here("img", "mainplot.png"), width = 900, height = 700)
(CO2_world | Top10)/ US_Indicators+
    plot_layout(widths = c(1, 2), heights = unit(c(4, 5), c('cm', 'null')))
dev.off()

quartz_off_screen 
                2